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Numerical Study on GRB-Jet Formation in CoUapsars 

Shigehiro Nagataki^'^, Rohta Takahashi^, Akira Mizuta^, Tomoya Takiwaki^ 

ABSTRACT 

Two-dimensional magnetohydrodynamic simulations are performed using the 
ZEUS-2D code to investigate the dynamics of a collapsar that generates a GRB 
jet, taking account of realistic equation of state, neutrino cooling and heating 
processes, magnetic fields, and gravitational force from the central black hole 
and self gravity. It is found that neutrino heating processes are not so efficient 
to launch a jet in this study. It is also found that a jet is launched mainly by 

fields that are amplified by the winding-up effect. However, since the ratio of 
total energy relative to the rest mass energy in the jet is not so high as several 
hundred, we conclude that the jets seen in this study are not be a GRB jet. This 
result suggests that general relativistic effects, which are not included in this 
study, will be important to generate a GRB jet. Also, the accretion disk with 
magnetic fields may still play an important role to launch a GRB jet, although 
a simulation for much longer physical time (~ 10 — 100 s) is required to confirm 
this effect. It is shown that considerable amount of ^^Ni is synthesized in the 
accretion disk. Thus there will be a possibility for the accretion disk to supply 
sufficient amount of ^^Ni required to explain the luminosity of a hypernova. Also, 
it is shown that neutron-rich matter due to electron captures with high entropy 
per baryon is ejected along the polar axis. Moreover, it is found that the electron 
fraction becomes larger than 0.5 around the polar axis near the black hole by 
Ue capture at the region. Thus there will be a possibility that r-process and 
r/p— process nucleosynthesis occur at these regions. Finally, much neutrons will 
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be ejected from the jet, which suggests that signals from the neutron decays may 
be observed as the delayed bump of the light curve of the afterglow or gamma- 
rays. 

Subject headings: gamma rays: bursts — accretion, accretion disks — black hole 
physics — MHD — supernovae: general — nucleosynthesis 



1. INTRODUCTION 



There has been growing evidence linking long gamma-ray bursts (GRBs; in this study, 
we consider only long GRBs, so we refer to long GRBs as GRBs hereafter for simplicity) 
to the death of massive stars. The host galaxies of GRBs are star -forming galaxies and 



the positions of GRBs appear to trace th e blue light of young stars (IVreeswijk et al. 2001 



Bloom et al. 20021 : iGorosabel et al. 20031 ). Also, 'bumps' observed in some afterglows can be 



naturally explained as contribution of brigh t supernovae (IBloom et al. 19991 : iReichart 1999 



Galama et al. 2000l : iGarnavich et al. 20031 ). Moreover, direct evidence of some GRBs ac- 
compani ed by supernovae ha.ve been reported such as the association of GRB 980425 with SN 



1998bw (IGalama et al. 19981: llwamoto et al. 19981) . that of GRB 030329 with S N 2003dh flHiorth et al. 200: 
Price et al. 20031: Istanek et al. 2003h . and that of GRB 060218 and SN 2006aj jMirabal et al. 20061 : 
Mazzah et al. 20061 ) . 



It should be noted that these supernovae (except for SN 2006aj) are categorized as a 
new type of supernovae with large kinetic energy (~ 1 0^^ ergs), nickel mass (~ O.SMq), and 
luminosity ( Iwamoto et al. 1998 : Woosley et al. 19991 ). so these supernovae are sometimes 
called hypernovae. The total explosion energy of the order of 10^^ erg is too important to be 
emphasized, because it is generally considered that a normal core-collapse supernova cannot 
cause such an energetic explosion. Thus another scenario has to be considered to explain 
the system of a GRB associated with a hypernova. One of the most promising scenarios is 
the collapsar scenario (jWoosley 19931 ). In the collapsar scenario, a black hole is formed as 
a result of gravitational collapse. Also, rotation of the progenitor plays an essential role. 
Due to the rotation, an accretion disk is formed around the equatorial plane. On the other 
hand, the matter around the rotation axis freely falls into the black hole. MacFadyen and 
Woosley (1999) pointed out that the jet-induced explosion along the rotation axis may occur 
due to the heating through neutrino anti-neutrino pair annihilation that are emitted from 
the accretion disk (see also Fryer and Meszaros 2000). 



It is true that the collapsar scenario is the breakthrough on the problem of the cen- 
tral engine of GRBs. However, there are many effects that have to be involved in order 
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to establish the scenario firmly. First of all, neutrino heating effects have to be investi- 
gated carefully by including microphysics of neutrino processes in the numerical simula- 
tions. It is true that MacFadyen and Woosley (1999) have done the numerical simulations 
of the collapsar, in which a jet is launched along the rotation axis, but detailed micro- 
physics of neutrino heating is not included in their simulations. Secondly, it was pointed 
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neutrino heating effects are not included in their works. Recently, Rockefeller et al (2006) 
presented 3-dimensional simulations of collapsars with smoothed particle hydrodynamics 
code. In their study, 3-fiavor flux-limited diffusion package is used to take into account neu- 
trino cooling and absorption of electron-type neutrinos, although neutrino anti-neutrino pair 
annihilation is not included. They have shown that alpha- viscosity drives energetic explosion 
through 3-dimensional instabilities and angular momentum transfer, although the jet is not 
launched and magnetic fields (source of the viscosity) are not included in their study. Thus 
it is not clear which effects are most important to launch a GRB jet, that is, what process 
is essential as the central engine of GRBs. 

Due to the motivation mentioned above, we have performed two-dimensional magne- 
tohydro dynamic simulations of collapsars with magnetic fields, rotation, and neutrino cool- 
ing/heating processes. In our simulations, the realistic equation of state (EOS) of Blinnikov 
et al. (1996) and effects of photo-disintegration of nuclei are also included. We investigated 
influence of magnetic flelds on the dynamics of collapsars by changing initial amplitude of 
the magnetic flelds. In section [21 models and numerical methods are explained. Results are 
shown in section [3l Discussions are described in HI Summary and conclusion are presented 
in section [51 



2. MODELS AND NUMERICAL METHODS 

Our models and numerical methods of simulations in this study are shown in this section. 
First we present equations of ideal MHD, then initial and boundary conditions are explained. 
Micro physics included in this study (equation of state (EOS), nuclear reactions, and neutrino 
processes) is also explained. 
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2.1. Magnetohydrodynamics 

We have done two-dimensional MHD simulations taking account of self-gravity and grav- 
itational potential of the central point mass. The calculated region corresponds to a quarter 
of the meridian plane under the assumption of axisymmetry and equatorial symmetry. The 
spherical mesh with 150(r)x 30(6') grid points is used for all the computations. The radial 
grid is nonuniform, extending from r =1.0x10^ to 1.0x10^'^ cm with finer grids near the 
center, while the polar grid is uniform. The minimum radial grid is set to be 3.0x10^ cm. 
We have confirmed that a free flow is smoothly solved with this grid resolution. 

The basic equations in the following form are finite differenced on the spherical coordi- 
nates: 
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where p, v, p, $p, $s e, L^, Lnuci, Ye, Yp, F„, Tp^n, and Tn^p, and are density, velocity, pres- 
sure, gravitational potential due to the central point mass (black hole), gravitational potential 
due to self gravity, internal energy density, heating/cooling rates due to neutrino processes, 
energy gain (loss) rate due to nuclear reaction, fraction of electron, proton, neutron, and 
reaction rate from proton to neutron (electron capture rate plus z/g capture on neutron) and 
from neutron to proton (positron capture plus Ue capture on proton), respectively. The La- 
grangian derivative is denoted as D/ Dt. The gravitational potenti al of the black hole is mod - 



ified to take into account some of the general relativistic effects fjPaczyhski fc Witta 1980l ). 
$p = —GM/{r — rs) where = 2GM/(? is the Schwartzshild radius. Self-gravity is obtained 
in ZEUS-2D code developed by Stone and Norman (1992a,b), by solving the matrix which 
results from finite-differencing the Poisson equation in two dimension. The ZEUS-2D code 
is also used to solve the MHD equations with second order accurate interpolation in space. 
Energy gain (loss) rate due to nuclear reaction and heating/cooling rates due to neutrino 
processes are described in subsections 12 . 3 . 2 | and 12 . 3 .31 Effects of a- viscosity and (anomalous) 
resistivity are not included in this study to avoid the uncertainty of the treatment of these 
effects. 
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2.2. Initial and Boundary Conditions 

We adopt the model E25 in Heger et al. (2000). This model corresponds to a star that 
has 25Mq initially with solar metallicity, but loses its mass and becomes to be 5A5Mq of a 
Wolf-Rayet star at the final stage. This model seems to be a good candidate as a progenitor 
of a GRB since losing their envelope will be suitable to be a Type Ic-like supernova and 
to make a baryon poor fireball. The mass of the iron core is l.GQM© in this model. Thus 
we assume that the iron core has collapsed and formed a black hole at the center. This 
treatment is same with Proga et al. (2003). The Schwartzshild radius of the black hole is 
5.0x10^ cm initially. 

We explain how the angular momentum is distributed initially. At first, we performed 
1-D simulation for the spherical collapse of the progenitor for 0.1 s when the inner most 
Si-layer falls to the inner most boundary (=10^ cm). After the spherical collapse, angular 
momentum was distributed so as to provide a constant ratio of 0.05 of centrifugal force to the 
component of gravitational force perpendicular to the rotation axis at all angles and radii, 
except where that prescription resulted in jiQ greater than a prescribed maximum value, 10. 
This treatment is similar to the one in MacFadyen and Woosley (1999). The total initial 
rotation energy is 2.44 x 10^^ erg that corresponds to 1.3 x 10~^ for initial ratio of the rotation 
energy to the gravitational energy (T/W). 

Configuration and amplitude of the magnetic fields in a progenitor prior to collapse 
are still uncertain. Thus in this study we choose a simple form for the initial configuration 
prior to collapse and the amplitude is changed parametrically. Initial configuration of the 
magnetic fields is chosen as follows: 

B(r) = (—) ^ (2 cos 0er + sin 0ee) for r > ro (7) 

3 \ r / 

2 ^ ^ 

= -BolcosOer — sm6eg) for r < ro. (8) 
3 

This configuration represents that the magnetic fields are uniform in a sphere (r < ro) , while 
dipole at outside of the sphere. We set ro to be the boundary between CO core/He layer (= 
3.6 x 10^ cm). Bq corresponds to the strength of the magnetic field in the sphere. We have 
chosen Bq to be 0, lO^G, lO^G, 10^°G, lO^^G, and IO^^q 

As for the boundary condition in the radial direction, we adopt the outfiow boundary 
condition for the inner and outer boundaries. That is, the flow from the central black hole 
is prohibited at the inner boundary and the inflow from the surface of the progenitor is 
prohibited at the outer boundary. Of course, the mass of the central black hole becomes 
larger due to the mass accretion from the inner boundary. As for the boundary condition 



- 6 - 



in the zenith angle direction, axis of symmetry boundary condition is adopted for the rota- 
tion axis, while the reflecting boundary condition is adopted for the equatorial plane. As 
for the magnetic fields, the equatorial symmetry boundary condition, in which the normal 
component is continuous and the tangential component is reflected, is adopted. 



2.3. Micro Physics 

2.3.1. Equation of State 

The equation of state (EOS) used in this study is the one developed by Blinnikov et 
al. (1996). This EOS contains an electron-positron gas with arbitrary degeneracy, which is 
in thermal equilibrium with blackbody radiation and ideal gas of nuclei. We used the mean 
atomic weight of nuclei to estimate the ideal gas contribution to the total pressure, although 
its contribution is negligible relative to those of electron-positron gas and thermal radiation 
in our simulations. 



2.3.2. Nuclear Reactions 

Although the ideal gas contribution of nuclei to the total pressure is negligible, ef- 
fects of energy gain/loss due to nuclear reactions are important. In this study, nuclear 
statistical equilibrium (NSE) was assumed for the region where T > 5 x 10^ [K] is satis- 
fied. This treatment is based on the assumption that the timescale to reach and maintain 
NSE is much shorter than the hydrodynamical time. Note that complete Si-burning oc- 
curs in explosive nucleosynt hesis of core-collapse supernovae for the region T > 5 x 10^ 



[K] (IThielemann et al. 19961 ). The hydrodynamical time in this study, ~ s (as shown in 
FigsO m and[TOl the nuclear reaction occurs at ~ (10^ — 10^)cm where the radial velocity 
is of the order of (10^ — 10®) cm s~^), is comparable to the explosive nucleosynthesis in core- 
collapse supernovae, so the assumption of NSE adopted in this study seems to be well. 5 
nuclei, n, p,^ He,^^ O, and ^^Ni were used to estimate the binding energy of ideal gas of nuclei 
in NSE for given (p, T, Ye). Ye is electron fraction that is obtained from the calculations of 
neutrino processes in section [2. 3. 3[ On the other hand, we assumed that no nuclear reaction 
occurs for the region where T < 5 x 10^ [K]. 
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2.3.3. Neutrino Processes 



Neutrino cooling processes due to pair capture on free nucleons, pair annihilation, and 
plasmon decay are included in this study. Since photoneutrino and b remsstrahlung pr ocesses 



are less important ones at 10^ < T < 10^^ [K] and p < 10^^ [g cm ^] (jitoh et al. 19891 ) where 



neutrino cooling effects are important in our calculations, we do not include these processes. 

As for the electron capture on free proton, we extend the formulation of Epstein and 
Pethick (1981) to arbitrary degeneracy. The energy loss rate per proton Qec [GeV s~^ 
proton"^ ] is 



Q 



EC 



X 



exp(F5(r/)/F4(r/)) 



where Gp = GpcosOc (in units of GeV^^) is the Fermi coupling constant [6^ is the Cabbibo 
angle), Cy = 1/2 + 2sin^6'iy is the vector coupling constant {9w is the Weinberg angle), 
Ca = 1/2 is the axial- vector coupling constant, is the Boltzmann constant, and Qy = 



\mr. 



m„ c 



is mass difference between neutron and proton, x and rj are defined as x 



(/ie — Qy)/{kBT) and rj = pe/iksT), where fie is the chemical potential of electrons. Fi{x) 
is the Fermi integral that is defined as 



Fix) 



1 + 



:dy. 



(10) 



We used the formulation derived by Takahashi et al. (1978) to solve the Fermi integral. 
Note that the leading term becomes the first term in Eq. [10] and obeys oc when elec- 
trons are not degenerated, although other terms dominate when electrons be come to be 
degenerated. We also include the Fermi blocking factor (IHerant et al. 19941 ) due to z/g 
assuming that the chemical potential of is zero and mean energy of the emitted neu- 
trinos by this process is kBTF^irj) / Fiiiji) (IRosswog &: Liebendorfer 20031 ). All we have to 
do is to change /ie to -/ie and to -Q^ to estimate the positron capture on free neu- 
tron. As for th e plasmon decays, we used the formulation of Itoh et al. (1989) with the 
blocking factor ( IHerant et al. 19941 ) assuming that the chemical potential of Ve is zero and 
mean energy of t he emitted neutrinos is k^T { 1 -|-7j^gg/2(l + 7ioss)} where 7ioss = 5.565 x 
10~^A/vr^/3 -|- rf ( Rosswog fc Liebendorfer 20031 ). We have also used the formul ation of Itoh 



et al. (1999) for pair annihilation with the blocking factor (IHerant et al. 19941 ). The mean 
energies of V f. and Pf. emitted by this process are kBTFiirf) / F^iirf) and kBTF4{—r])/F3{—r]), 
respectively (ICooperstein et al. 19861 ). 
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In our simulations, density at the inner most region reaches as high as 10^^ g cm~^ 
where nuclei are almost photo-disintegrated into nucleons. In this region, as discussed in 
section 14.11 optically thin limit breaks down. In ord er to take into acc ount such a high 



density region, neutrino leakage scheme is introduced (iKotake et al. 20031 ). In this scheme, 
we calculate typical diffusion timescale for neutrinos to escape from the system. Since 
chemical composition is dominated by nucleons at such a high density region, the opacity is 
mainly determined by scattering and absorption by nucleons for v^, and z/g and scattering by 
nucleons for i/^, i/^, z/,-, and z/^. The mean free path can be written as 

K} = — o-sc(e^J + — o-ab(e,.J for (H) 
rriu ruu 

K} = — crsc(epj + — cTab(epJ for z/e (12) 
rriu rriu 

Aj7^ = -^Osci^vx) for other flavors (13) 

rriu 

where m„, a^c, CTab, Yp, Yn, e^^, ep^, and e^^ are the atomic mass unit, scattering cross section 
by nucleon, absorption cross section by nucleon, mass fraction of proton and neutron, and 
mean energies of u^, z/e, and u^. We take the cross sections for these interactions from Bethe 
(1990). The diffusion timescale is given by 

where we take AR = (rArA6')^/^ as the typical diffusion length. Ar and rA9 are grid sizes 
in the radial and polar directions at that point. Neutrino emissivity [GeV cm~'^ s~^] at the 
point is calculated as 

/^'^ = /,xminfl.O,— ) (15) 

where f^, and dt are neutrino emissivity obtained by the optically thin limit and time step 
of the simulations. On the other hand, when we estimate the neutrino heating rate, we take 
optically thin limit for all the regions for simplicity. This corresponds to that we estimate 
the upper limit for the neutrino heating effects. Even this limit, as shown later, jet like 
explosion does not occur by the neutrino heating (as shown in Fig. [1]). Also, this optically 
thin limit will be justified since neutrino cooling rate dominates heating rate, as shown in 
Fig. El Detailed discussion on this point is presented in section 14. 1[ 

Neutrino heating process due to neutrino pair annihilation and i^e and z/g captures on 
free nucleons with blocking factors of electrons and positrons are included in this study. As 
for the neutrino pair annihilation process, the formulation of Goodman et al. (1987) (with 
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blocking factors) is adopted. The and captures on free nucleons are inverse processes 
of electron/positron captures. The calculation of neutrino heating is the most expensive in 
the simulation of this study. Thus, to save the CPU time, the neutrino heating processes 
are calculated only within the limited regions (r < 10^ cm). Moreover, we adopt some 
criterions as follows in order to save CPU time. We have to determine the energy deposition 
regions and emission regions to estimate the neutrino heating rate. As for the neutrino pair 
annihilation process. We adopt no criterion for the energy deposition regions other than 
r < 10^ cm, while the emission regions satisfy the criterion of T > 3 x 10^ K. As for the 
neutrino capture process, we adopt the criterions as follows: The absorption regions should 
satisfy the criterion of p > 10^ g cm~^, while the emission regions satisfy the criterion of (i) 
p > 10^ g cm~^ and (ii) T > 10^ K. Also, these heating rate had to be updated every 100 
time steps to save CPU time. Of course this treatment has to be improved in the future. 
However, this treatment seems to be justified as follows: The total time step is of the order 
of 10^ steps. The final physical time is of the order of seconds, so the typical time step is 
~ 10^^ s. This means that heating rate is updated every ~ 10~^ s, which will be shorter 
than the typical dynamical timescale. This point is discussed in detail in section 14.11 

The heating process of neutrino pair annihilation are calculated as follows. Let points 
A and B be in the emission regions and a point C be in the energy deposition regions. Let 
c/a, ds-, and dc be distances between A and B, B and C, and C and A, respectively. The 
angle is defined as the angle ACB. The energy deposition rate at C, Q{C) [GeV cm"^ s~^], 
is given formally by integrating the emission regions as 

Q(C) = ^ / / ^4^4 ^^ "^^2"^^' / rfe.^^epe,ep(e, + ep)F,Fp(l-/-)(l-/+), (16) 

where is (1—4 sin^ + 8 sin^ Qw)l (6??) for t/^,z/^7 and v^-v^ annihilation and (1 + 4 sin^ ^ + 
8 sin'' ^vt/)/(67r) for v^aj^ annihilation (iGoodman et al. 19871 ). xa and xb are locations of the 



points A and B, e^/^p are energies of (anti-)neutrinos, F^^y are energy spectrum [cm~^ s~^ 
GeV"^] of (anti-)neutrinos, and f"^ are Fermi blocking factors of electrons/positrons at point 
C. In this formulation, Gp is written in units of cm GeV^^. It is easily seen that Eq. ( |T6l) 
is eight dimensional integral and it takes too much time to accomplish this integral. Thus 
we have approximated that each neutrino energy spectrum due to each emission process 
is a monotonic one, that is, we assumed only neutrinos with average energy are emitted 
when Eq. (fT6l) is carried out. Since the energy loss rate [GeV cm~^ s~^] is obtained from 
the formulations of neutrino cooling process mentioned above, number of emitted neutrinos 
for each process is calculated by dividing the average energy of emitted neutrinos for the 
process. Although the calculated region is a quarter of the meridian plane, the neutrino 
flux from the south hemisphere is included assuming the symmetry relative to the equatorial 
plane. The same treatment with neutrino pair annihilation process is adopted to estimate the 
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heating rate due to Ve and capture processes, that is, average energy of (anti-)electron-type 
neutrinos is used. 

Neutrinos are emitted isotropically in the fluid-rest frame, so strictly spea king, neutrinos 



are emitted unisotropically in the coordinate system due to the beaming effect (IRybicki and Lightman 19791 ! 

In fact, the angular frequency is found to become as large as 10'^ s~^ at the inner most region 

(Figs. [3] and [TOl) and the rotation velocity, reaches ~ 10^° cm s~^ at most. Thus the 

beaming effect may be important although we did not take into account the effect in this 

study. 

Finally the models considered in this study are summarized in Table [TJ The digit in the 
name of each model represents the power index of Bq. 



3. RESULTS 

We present results of numerical simulations in this section. First of all, we show in Fig[T] 
the density contour of the central region of the progenitor (r < lO^cm) with velocity fields 
for model and model 9. In the top panels, density contour with velocity fields for model 
at t = 2.1 s (left panel) and t = 2.2 s (right panel) is shown, while in the bottom panels 
the ones for model 9 is shown. The color represents the density (g cm~^) in logarithmic 
scale (10^ - 10^2). Vertical axis and horizontal axis represent polar axis (= rotation axis) 
and equatorial plane, respectively. You can easily find that a jet is launched at t = 2.1 s 
for model 9, while no jet occurs for model 0. In the following subsections, we explain the 
dynamics of these models in detail. 



3.1. Dynamics without Magnetic Fields 

In this subsection, we mainly explain the dynamics of model 0, then the one of model 9 is 
explained in the next subsection. First, we show the accreted energy and mass accretion rate 
as a function of time in Fig. [2]^a). Solid, short-dashed, and dotted lines represent accreted 
total energy, kinetic energy, and thermal energy ( x 10^° erg) into the black hole. Long-dashed 
line represents the mass accretion rate {Mq s~^) as a function of time. From this figure, we 
can understand the following points: (i) the accreted energy amounts to of the order of 10^^ 
erg, which is comparable to the explosion energy of a hypernova. Thus there is a possibility 
for the models in this study to explain the large explosion energy of a hypernova as long as 
the released gravitational energy is efficiently converted to the explosion energy, (ii) kinetic 
energy dominates thermal energy at the inner boundary. This is because almost all thermal 



- 11 - 



energy is extracted in the form of neutrinos, which is seen in Fig. [6]^a). (iii) mass accretion 
rate drops almost monotonically from ~ 10^^ Mq to ~ IO^^Mq s^^, which is discussed 
in sections 13.21 and 1 4.21 

Next, we show in Fig. [3] profiles of physical quanta of the accretion disk around the 
equatorial plane. Profiles of density, absolute value of radial velocity, angular frequency, 
temperature, density scale height, and specific angular momentum in the accretion disk are 
shown for model at t = 2.2 s. From the figure, we can understand the following points: (i) 
the density reaches as high as 10^^ g cm~^ around the central region, (ii) the inflow velocity 
{vr) becomes as low as 10^ cm s~^ at the central region, by which the mass accretion rate 
becomes as low as 10~^ Mq s~^, as shown in Fig. [2](a) (Note that the mass accretion mainly 
comes from the region between the rotational axis and the disk. See also Igumenshchev and 
Abramowicz 2000; Proga and Begelman 2003). 

We show in Fig. H] profiles of mass fraction for nuclear elements at the equatorial plane 
for model at t = 2.2 s. Dot-dashed, dotted, short-dashed, long-dashed, and solid lines 
represent mass fraction of n, p, ^He, ^^O, and ^^Ni, respectively. From this figure, we can 
easily see that oxygen is photo-dissociated into helium at 2 x lO^cm, while helium is photo- 
dissociated into nucleons at 3 x lO^cm. It is also noted, some ^^Ni is seen at r < 3 x 10^ cm 
(solid line), which may explain the luminosity of a hypernova as long as it is ejected. This 
point is discussed in section 14. 3[ The discontinuity of temperature at 3 x 10^ cm in Fig. [3] 
will come from the cooling effect due to the photo-disintegration of helium into nucleons. 

From now on we show the results on neutrino processes. In Fig. [5l we show contour 
of neutrino cooling rate with velocity fields for model at t = 2.2 s. The color represents 
the emissivity of neutrinos (erg cm~'^ s~^) in logarithmic scale (10^° — 10'^^). This emissivity 
of neutrinos is almost explained by pair-captures on free nucleons, as shown in Fig. O We 
can easily see that emissivity of neutrinos is high at the region where the accretion disk is 
formed. 

We show the results on neutrino cooling for every neutrino process. We show in Fig. [6t^a) 
cumulative energy ( x 10^° erg) of emitted neutrinos for each process as a function of time for 
model 0. Dot-dashed, long-dashed, solid, short-dashed, and dotted lines represent plasmon 
decay, electron-positron pair annihilation, positron capture, electron capture, and summation 
of all processes. It is clearly seen that almost all emitted energy comes from pair captures on 
free nucleons. Also, the emitted energy amounts to of the order of 10^^ erg (strictly speaking, 
3.44 X lO^^erg), which is comparable to the accreted kinetic energy and much higher than the 
accreted thermal energy (see Fig. [21^a)). Thus we consider that almost all thermal energy, 
which was comparable to the kinetic energy in the accretion disk, is extracted by neutrino 
emission. Note that time evolution of electron fraction is mainly determined by positron and 



- 12 - 



electron capture processes in the accretion disk. 

In Fig. [71 we show contour of neutrino heating rate with velocity fields for model at t 
= 2.2 s. The color represents the energy deposition rate (erg cm~'^ s~^) in logarithmic scale 
(10^° — 10^^). Left panel shows the energy deposition rate due to z/g and Ve captures on free 
nucleons, while right panel shows the energy deposition rate due to v and v pair annihilation. 
In the pair annihilation, contributions from three fiavors are taken into account. It is, of 
course, that contour of energy deposition rate due to z/g and v^. captures traces the number 
density of free nucleons, so the energy deposition rate is high around the equatorial plane 
where the accretion disk is formed. On the other hand, energy deposition rate due to v and v 
pair annihilation occurs everywhere, including the region around the polar axis. This feature 
will be good to launch a jet along the polar BjXIS, clS pointed by MacFadyen and Woosley 
(1999). However, this heating effect is too low to launch a jet in this study (see Fig. [T]). 

In Fig. [U^a), we show neutrino luminosity (solid line), energy deposition rate due to v 
and V pair annihilation (dotted line), and energy deposition rate due to z/g and v^. captures 
on free nucleons (dashed line) as a function of time for model 0. It is clearly seen that energy 
deposition rate is much smaller than the neutrino luminosity, which supports our assumption 
that the system is almost optically thin to neutrinos. Also, we can see that v^. and z/g captures 
on free nucleons dominates v and v pair annihilation process as the heating process. It is 
also noted that neutrino luminosity and energy deposition rate decreases along with time, 
which refiects that the mass accretion rate also decreases along with time (Fig. [2|^a)). 

Finally, we show in Fig. [9](a) the integrated deposited energy (xlO^° erg ) due to v 
and V pair annihilation (solid line) and z/g and z/g captures on free nucleons (dashed line) 
as a function of time. It is confirmed that z/g and Ve captures on free nucleons dominate 
V and V pair annihilation process as the heating process. However, as shown in Figs. 
and [71 z/g and v^. captures on free nucleons occur mainly in the accretion disk, where neutrino 
cooling effect dominates neutrino heating effect. Thus z/g and z/g captures on free nucleons 
is not considered to work for launching a jet. As for the v and v pair annihilation process, 
although this process deposit energy everywhere, including the region around the polar axis, 
the deposited energy amounts to only of the order of 10^^ erg. This is 10~^ times smaller 
than the explosion energy of a hypernova. In fact, as shown in Fig. [H the jet is not launched 
in model 0. Thus we conclude that the efficiency of neutrino heating is too low to launch a 
jet in this study. 
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3.2. Dynamics with Magnetic Fields 

In this subsection, we explain the dynamics of model 9 as an example of a collapsar 
with magnetic fields. Dependence of dynamics on the initial amplitude of magnetic fields is 
shown in section [3 .3 .11 

First, the accreted energy and mass accretion rate as a function of time are shown in 
Fig. [2]^b). The meaning of each line is same with Fig. [2](a), although dot-dashed line, which 
is not in Fig. [2] (a), represents accreted electro-magnetic energy (xlO^° erg) into the black 
hole as a function of time. The reason why the final time of the simulation for model 9 is 
2.23 s is that the amplitude of the magnetic field (in particular, B^) becomes so high at the 
inner most region that the Alfven crossing time at the region make the time step extremely 
small. Also, it seems that the mass accretion rate does not decrease so much in model 9. 
Rather, it seems to keep ~ O.O5M0 s~^. We guess this is because magnetic fields play a role 
to transfer the angular momentum from inner region to outer region, which makes matter 
fall into the black hole more efficiently. This point is also discussed with Figs. 151— [T^ 

Next, we show in Fig. [10] profiles of physical quanta of the accretion disk around the 
equatorial plane for model 9 at t = 2.2 s. When we compare these profiles with the ones in 
Fig. [31 we can see that the radial velocity is higher in model 9 at small radius. This may 
refiect that the mass accretion rate is higher in model 9 than model (as stated in section 
3.1, the mass accretion mainly comes from the region between the rotational axis and the 
torus. So we have to note that this feature does not explain the mass accretion rate directly, 
see also Fujimoto et al. 2006). As for the other profiles, there seems no significant difference 
between the two models. 

In Fig. [m we show profiles of amplitude of magnetic fields as a function of radius on 
the equatorial plane for model 9 at t = 2.2 s. Dotted, solid, and dashed lines represent 
the amplitude of B^, Bg, and B^f,, respectively. It is clearly seen that i?^ dominates within 
r = lO^cm. Thus we can conclude that magnetic pressure from B^ drives the jet along the 
rotation axis (see Fig.[T]). This point is also discussed with Figs. [T31— [T31 Also, these magnetic 
fields may play a rol e to transfer the angular momentum. Since t he viscosity parameter a 



can be estimated as (iBalbus fc Hawley 19981 : lAkiyama et al. 20031 ) 



BrB^ 

" - (17) 

we plot in Fig. [11] the estimated viscosity parameter (xlO^^'^) (long-dashed line). This figure 
suggests that a-viscosity coming from the magnetic fields may play a role to transfer the 
angular momentum at inner most region effectively. However, the angular momentum cannot 
be transfered to infinity along the radial direction. This is confirmed by the Alfven Mach 
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number (xlO~^°) in the radial direction (= Vr/vr^A] dot-dashed hne in Fig. [TT]) . At only 
the inner most region, the flow becomes marginally sub-Alfvenic where the viscous force 
due to magnetic stress can bring the angular momentum outward. Thus we consider that 
the outflow (including the jet) in the polar direction (see Fig. [T] bottom) should bring the 
angular momentum from the inner most region. 

As additional information, we found that the velocity of the slow magnetiosonic wave is 
almost same with the Alfven velocity. On the other hand, we found that the fluid is subsonic 
against the fast magnetiosonic wave in the simulated region. 

We show in Fig.[T2]the proflles of speciflc angular momentum (cm^ s~^) on the equatorial 
plane for model (dotted line) and model 9 (short-dashed line) at t = 2.2 s. The proflles of 
angular momentum density (g cm~^ s~^) for model (long-dashed line) and model 9 (solid 
line) are also shown in the flgure. We conflrmed that the speciflc angular momentum is 
not so different between model and model 9, but it is found that the angular momentum 
density is lower in model 9 compared with model at the inner region. This feature might 
reflect that the matter falls into the black hole efficiently in model 9 at small radius (see 
Fig. [TT]) . This picture seems to be consistent with the almost constant accretion rate in 
model 9 (Figs. [2](b)). 

Cumulative energy (xlO^° erg) of emitted neutrinos for each process as a function of 
time for model 9 is shown in Fig. E^b), neutrino luminosity, energy deposition rate due to v 
and v pair annihilation, and energy deposition rate due to z/g and z/g captures on free nucleons 
as a function of time for model 9 are shown in Fig.[8](b), and integrated deposited energy due 
to V and v pair annihilation and z/g and Ve captures on free nucleons as a function of time for 
model 9 are shown in Fig. [9](b) The meaning of each line in each flgure is same with the one 
used for model 0. We can derive a similar conclusion for the role of neutrino heating effect 
in model 9 as in model 0. As for the v^, and z/g captures on free nucleons, the cumulative 
deposited energy becomes as high as 10^^ ergs that is comparable to the explosion energy of a 
hypernova. However, this heating process occurs mainly in the accretion disk, where neutrino 
cooling effect dominates neutrino heating effect. Thus z/g and z/g captures on free nucleons is 
not considered to work for launching a jet. As for the v and v pair annihilation process, the 
deposited energy amounts to no more than 10^° erg. Thus we conclude that the neutrino 
heating effect in model 9 is too inefficient to launch a GRB jet and cause a hypernova. It has 
to be noted that the energy deposition rate due to pair captures on free nucleon sometimes 
becomes larger than the neutrino luminosity in Fig. [8](b). This means that the optically 
thin limit breaks down at that time. This point is also discussed in section 14.11 We also 
found that the neutrino luminosity, energy deposition rate, and integrated deposited energy 
seem to be higher in model 9 than in model 0. We consider that this feature comes from 
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the high mass accretion rate (high rate of release of gravitational energy) caused by angular 
momentum transfer due to magnetic fields. 

We show in Fig. [T3]time evolution of energy of magnetic fields (xlO^° erg) in the whole 
calculated region for the case of model 9. Dotted, solid, and dashed lines represent energy in 
the form of Br, Be, and B^, respectively. We can see that B^ grows rapidly and dominates 
other components, although B^ is much smaller than other components at first (note that 
1-D simulation of the spherical collapse of the progenitor is done until t =0.1 s (see section 
2.2), so B^ is set to be until t = 0.1 s). i?^ is amplified by winding effect and launches 
a jet along the polar axis (see Fig. [1]). The total energy of 5^ at the final stage of the 
simulation in this study is of the order of 10^° erg, which is much smaller than the explosion 
energy of a hypernova. However, the total energy of B^ is increasing almost monotonically. 
As mentioned above, the final time of the simulation for model 9 is determined due to the 
reason that the Alfven crossing time becomes extremely small at the inner most region. 
When we can overcome this problem and simulate the model 9 further, there is a possibility 
that the energy of B^ becomes comparable to the explosion energy of a hypernova. The 
evolution of Br and Bg is similar to the one by Magnetic Rotational Instability (MRI), 
which is discussed in section 14.21 

We show in Fig. [H] contour of plasma beta (=Pgas+radiation/Pmag) with magnetic fields 
{Br and Bg) for model 9 at t = 2.2 s. In Pgas+radiation, degenerated pressure of electrons 
is included. The color represents the plasma beta in logarithmic scale (10~^ — 10®). The 
minimum value of plasma beta in this region is 0.193. We can see clearly that the beta value 
is low in the jet region, from which we can understand the jet is launched by the magnetic 
fields (in particular, B^). 

Finally, we have performed these simulations presented in sections 13.11 and 13.21 without 
calculating self-gravity, and found that the dynamics of the coUapsar is not changed so much. 

3.3. Toward Discussions 

Before we discuss the results mentioned above, we present two more results in this 
section. One is the dependence of dynamics on initial amplitude of magnetic fields, the 
other is the explosive nucleosynthesis. 
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3.3.1. Dependence on Initial Amplitude of Magnetic Fields and on Resolution of Grids 

Here we show the dependence of dynamics on initial amphtude of magnetic fields. In 
Fig. [151 evolutions of total energy of magnetic fields (xlO^° erg) for model 8 (top- left panel), 
model 10 (top-right panel), model 11 (bottom-left panel), and model 12 (bottom right panel) 
are shown. The final time of the simulation is determined, like in model 9, by the reason that 
the Alfven crossing time becomes extremely small at the inner most region. In all models, 
the energy of Bfp amounts to 10^^ — 10^° erg, which means that the amplitude of becomes 
as strong as lO^^G. This can be understood by simple calculations: the typical radius of the 
inner most region is of the order of 10^ cm, so the volume of the region times the energy 
density of magnetic fields becomes Eb = 8.3 x 10^^{B/2 x 10^^G)^(r/5 x lO^cm)^ erg. 

The total energy of Br and Bg components does not change so much with time. In 
particular, in models 10 — 12, their total energy hardly changes. When we compare Br with 
Bg, Br seems to be more unstable than Bg (see also Fig. [T^ . This is similar to the results 
of the local simulations of MRI, which is discussed in section 14.21 

We summarize the features of the launched jet at the final stage of simulations for every 
model in Table [2l We define the jet as follows: (i) the region within some opening angle {6j) 
(ii) the region where total energy (i.e. summation of kinetic, thermal, electro-magnetic, and 
gravitational energies) is positive at the final stage of the simulations, (iii) the amplitude 
of the velocity is larger than 5 x 10^ cm s~^. In Table |2l total energy, mass, terminal 
bulk Lorentz factor, and ratio of the magnetic energy relative to total energy of the jet is 
shown assuming that 6j = 5°, 10°, and 15° for each model. Terminal bulk Lorentz factor is 
estimated by assuming that total energy goes into kinetic energy during expansion. 

We can see that terminal bulk Lorentz factor in every model is much smaller than the 
required value (~ 300) for GRB jets. Also, the energy of the jet in every model is much 
smaller than the typical energy of a GRB. Thus we conclude that the jets seen in this study 
will not be GRB jets. 

Also, we show in left panel of Fig. [16] the evolution of total energy of magnetic fields 
for the case of model 9 (same with Fig. [T3|) . but for 150(r)x5(^^) grid points (left panel) 
and 150(r) x20(6') grid points (right panel) to show the dependence of results on the grid 
resolution. It is found that the growth rate of Br and Bg components depends on the grid 
resolution. On the other hand the growth rate of B^ does not depend on the grid resolution 
(that is, the energy in B^ fields becomes 10'^^ — 10^° ergs in (1-2) sec). The saturation level of 
Br and Bg does not depend on the grid resolution so much. This is discussed in section 14.21 
We have also done a simulation of model 9 with finer resolution (300(r) x60(6')), although 
neutrino anti-neutrino pair annihilation is not included to save CPU time. The simulation 
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region is set to be (10^ cm < r < 10^ cm, < 6 < 90). The minimum radial grid is set to be 
3x 10^ cm. We have found that the dynamics of collapsars is hardly changed. We have found 
that an accretion disk is formed around the black hole and a jet is launched at t = 1.98 sec. 
Moreover, we have found that the evolution of magnetic fields is hardly changed, which is 
shown in right panel of Fig. [161 This result means that a jet is driven by magnetic fields and 
the standard resolution of our study (150(r) x30(6')) is not so bad. 

3.3.2. Nucleosynthesis 

In this section, we present results on explosive nucleosynthesis. We show in Table [3] 
the mass of ^^Ni in the regions where total energy is positive {M^^ (^0)) ^ind the total 
mass of ^^Ni in the whole simulated region (Mj^°* (Mq)) for each model. These amounts are 
estimated at the final stage of the simulations. The ejected mass of ^^Ni, M^^, is too little to 
explain the luminosity of a hypernova, although considerable amount of ^^Ni is synthesized 
in the accretion disk (M^°* (Mq)). We found that the ejected mass is mainly composed of 
n, p, and He. This is supported by Fig. [T71 where entropy per baryon at t = 2.2 s for model 
9 is shown in units of ki,. In the jet region, the entropy per baryon is remarkably high, so it 
is natural that light elements dominate in the jet region. 

Finally, we show in Fig. [TS] contour of electron fraction (Yg) with velocity fields for 
model at t = 2.2 s (left panel) and model 9 at t = 2.2 s (right panel). The color represents 
the electron fraction on a linear scale (0.1-0.540). We can easily see that Ye becomes low in 
the accretion disk. This is because electrons are degenerated and electron capture dominates 
positron capture at this region. Also, we can see that mass element with low Yg is ejected in 
model 9. Also, mass element with high (highest value is 0.522) is also ejected in model 9 
from the inside of the low Ye jet along the polar axis near the black hole. 

Since entropy per baryon is very high in the jet region (Fig. [TT]) . these mass elements 
may cause r— process and/or r/p— process nucleosynthesis. Here we have to comment on the 
electron fraction at the high density region. The value of the electron fraction is solved of 
the order of 0.1 in the accretion disk. This value is obtained assuming that the chemical 
potential of electron-type neutrinos is zero. If careful treatment of neutrino transfer is done, 
the chemical potential of electron-type neutrinos may prohibit the electron fraction from 
being as low as 0.1. This point is discussed in sections Wl] and 
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4. DISCUSSIONS 

In this section, we discuss our numerical results and prospect for future works. We 
discuss neutrino physics, effects of magnetic fields, nucleosynthesis, general relativistic effects, 
initial conditions, and prospect for improvements of our numerical code. 



4.1. Neutrino Physics 

In this study, we found that deposited energy due to neutrino pair annihilation are too 
small to explain the explosion energy of a hypernova and a GRB (Fig. [9]). Even though 
the deposited energy by electron-type neutrino capture on free nucleons can be comparable 
to the explosion energy of a hypernova in model 9 (Fig. [H]), the deposition region is the 
high density region (Fig. [7]) where cooling effect dominates the heating effect (Fig. In 
particular, no jet was found in the numerical simulations of model (Fig. [1]). The energetics 
of this system can be understood from Figs. [2], [6l and [91 The released gravitational energy 
by collapse is the energy source of accreted energy and neutrinos. From Figs. [2] and [6], we 
can understand that the kinetic energy and thermal energy share the released gravitational 
energy almost equally, then almost thermal energy was extracted in the form of neutrinos. 
The total energies of accreted energy and neutrinos are of the order of 10^^ erg, then the 
emitted neutrino energy was deposited into matter through weak interactions. Its efficiency 
is less than 1% for neutrino pair annihilation, as can be seen in Figs. [6] and [91 As for the 
efficiency of electron-type neutrino capture it amounts to ~ 10 — 20%. This means that 
the inner most region of the accretion disk becomes to optically thick against neutrinos and 
~ 10 — 20% of neutrinos are absorbed. 

The deposited energy by neutrino pair annihilation is of the order of (10^^ — 10^° erg, 
which is much smaller than the explosion energy of a hypernova and a GRB. In order to 
enhance the deposited energy by neutrino pair annihilation, there will be two ways. One 
is to enhance the released gravitational energy, the other is to enhance the efficiency of 
energy deposition. The former corresponds to enhance the mass accretion rate, which will 
be realized if effective angular momentum transfer is realized. From Figs. [21 [HI and [T21 it 
was inferred that magnetic fields seem to work efficiently so that high mass accretion rate 
is realized. Of course, the mass accretion rate also depends on the distribution of initial 
angular momentum. We should investigate these effects further in the future. As for the 
efficiency of energy deposition, it will be enhanced when the general relativistic effects are 
taken into account. This is because neutrinos are trapped around the black hole, so that the 
possibilities of neutrino pair annihilation and neutrino capture become enhanced. This effect 
is investigated in detail by using a steady solution of an accretion disk in the forth-coming 
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paper. Of course, we are planning to include this effect in our numerical code in the future. 

Although we believe that our conclusion on the energetics mentioned above will be 
unchanged, we have to improve our treatment on the neutrino heating for further study. In 
this study, we took the optical thin limit to estimate the neutrino heating rate. This will 
be justified by Figs. [6]and[9l However, for further study, we have to investigate the cases in 
which mass accretion rate is higher than in this study to achieve energetic explosion enough 
to explain the explosion energies of a hypernova and a GRB. In fact, we consider that the 
optically thin limit breaks down even the models in this study at the highest density region. 
This is estimated as follows: The cross section of z/g and Ve captures on free nucleons is given 
by 0" ~ ao^ey/meC^Y where cxo = 1.76 x 10~^^ cm^. Since the highest density in the accretion 
disk is of the order of 10^^ g cm~^ at r ~ 10^ cm (the scale hight is also of the order of 
lO^cm) and typical energy of neutrinos are of the order of 10 MeV (see Figs. 151 and [TU]) . the 
optical depth at this region is r = a{p/mp)L ~ 4.2(ei./10MeV)^(p/10^^g cm"^)(L/10''cm). 
Thus, at the highest density region, the optically thin limit must break down. This picture is 
also confirmed in Fig. [8](b). In Fig. Mjo), as stated in section [3^ the energy deposition rate 
due to pair captures on free nucleon sometimes becomes larger than the neutrino luminosity. 
This reflects that the optically thin limit breaks down at that time. Although we believe that 
our conclusion on the energetics will not be changed so much, we are planning to develop 
the careful neutrino transfer code that includes emissions, absorptions, and scattering of 
neutrinos for further study. We also note four points that have to be improved for the 
treatments of neutrino heating. One is that we did not take into account the light-crossing 
time of the system and assumed that the system is almost steady during the light-crossing 
time when we estimate the neutrino heating rate. From Figs. [5] and [TJ the neutrino cooling 
and heating occur efficiently within several times of lO^cm. Thus the typical light crossing 
time will be of the order of 1ms. For comparison, the rotation period at the inner most region 
is ~ 6.3 X 10~^ s. Since the system forms an accretion disk and the viscosity parameter a 
is 0.1 at most (Fig. [TT]), the system will be treated steady at least ten times of the rotation 
period, 6.3 x 10~'^ s. Thus the treatment to neglect the light crossing time will be fairly 
justified. Second is that we update the neutrino heating rate every 100 timesteps to save 
the CPU time. The inner most radius is set to be 10^ cm, so the typical time step is 
estimated to be 10^ x A6/c s, where c is the speed of light and A9 = 7r/60. Thus 100 
time steps corresponds to 1.74 ~ 10~^ s, which will be comparab le to free-fall timescale 



{ts = 1/V2AttGp ~ 4.5 x IQ-^IO^^ g cm-y py/'^ s flWooslev 19861 ) ) and rotation period. 
Thus we believe that this treatment will be fairly justified. Third, we have approximated 
that each neutrino energy spectrum due to each emission process is monotonic, that is, 
we assumed only neutrinos with average energy are emitted when Eq. f ll6p is carried out. 
However, the cross section of neutrino pair annihilation is proportional to square of the total 
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energy in the center of mass, and that for electron-type neutrino absorption on free nucleons 
is proportional to square of neutrino's energy. Thus the contribution of neutrinos with 
high energy will enhance the efficiency of neutrino heating. These points will be improved 
when we can include a careful neutrino transfer code in future. The last point is related 
with the nuclear reactions. In this study, the NSE was assumed for the region where T > 
5 X 10^ [K]. Thus, the reactions to maintain NSE occur suddenly when the temperature 
becomes so high as to satisfy the criterion. However, in reality, NSE might break down at 
low density region, where cooling effect due to photo-disintegration will be not so strong 
as in this simulation. It should be also noted that the photo-dissociation from He into 
nucleons is strong cooling effect and absorb thermal energy when this reaction is switched 
on. Thus the thermal energy suddenly absorbed by nuclear reactions. That is seen as the 
discontinuities of temperature in Figs.[3]and[Tni Although we believe these discontinuities do 
not change our conclusion on the energetics mentioned above (because much more neutrinos 
come from the inner region; Fig. [5]), the profile of temperature will be solved smoothly when 
we use a nuclear reaction network instead of using the NS E relation. Since the emissivity 



of neutrinos depends very sensitively on the temperature (IBethe 1990l : iHerant et al. 1992 



Lee fc Ramirez- Ruiz 20061 ). estimation of temperature should be treated carefully. We are 
planning to check the dependence of temperature on the nuclear reaction network and several 
EOS in the forth coming paper. 

Finally, we discuss the detectability of neutrinos from collapsars. Since the event rate is 
much smaller than the normal core-collapse supernova, the chance probability to detect neu- 
trino signals from a coUapsar will be very small. However, if it occurs nearby our galaxy, the 
neutrino signal from a collapsar will be distinguished from normal core-collapse supernova. 
As for the normal core-collapse supernovae, the time evolution of the luminosity of neutrino 
of each fiavor is determined firmly by the binding energy of a neutron star and opacity of 
neutron star against neutrinos. On the other hand, in the case of a collapsar, the time 
evolution of neutrino luminosity will depend on the time evolution of mass accretion rate, 
which in turn should depend on the initial distribution of angular momentum and magnetic 
fields. Thus there should be much varieties of time evolution for the luminosity of neutrinos 
in the case of collapsars. Also, in the case of collapsars, the dominant process to generate 
neutrinos is pair captures on free nucleons (see Fig. [H]), so in the case of a collapsar, the 
electron-type neutrinos will be much more produced compared with other fiavors. This is in 
contrast with the normal core-collapse supernovae(e.g. Buras et al. 2006, and see references 
therein). Of course, we have to take vacuum and matter oscillation effects into account to 
estimate the spectrum of neutrinos from a collapsar precisely. In particular, in the case of a 
collapsar, the density distribution is far from spherically symmetric, so we have to be careful 
about viewing angle to estimate the matter oscillation effect. It is true that the event rate 
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of collapsars is smaller than normal core collapse supernovae, but the released gravitational 
energy can be larger if considera ble amount of mass of the progenitor falls into the central 
black hole (iNagataki et al. 20021 ). Thus we consider that there will be also a possibility to 
detect a neutrino background from collapsars. 



4.2. Effects of Magnetic Fields 

We have seen that the mass accretion rate seems to be enhanced in model 9, compared 
with model (Fig. [2]), which enhances the luminosity of neutrinos (Fig. [6]) and energy depo- 
sition rate due to weak interactions (Fig. [H]). This seems to be because magnetic viscosity is 
effective at inner most region (Fig. [11] and Fig. [T2|) and multi-dimensional outflow (Fig. [1]) 
carries angular momentum outward. Thus amplification of magnetic fields is important not 
only for launching a jet by magnetic pressure, but also for enhancing mass accretion rate and 
energy deposition rate through weak interactions. In our study, we assume the axisymmetry 
of the system to save CPU time. In this case, the field built up by the effec t of magnetoro - 



tational instability (MRI) decays due to Cowling's anti-dynamo theorem (jShercliff 19651 ). 
However, the plasma beta becomes lower than unity in the jet region (Fig. [T^ . which is 
embodied by the amplification of fields (Figs. [T3] and [T5|) . Thus we consider that B^ 
field is not amplified by MRI effects, but by winding-up of poloidal fields due to differential 
rotation. The typical timescale of winding-up at the inner most region will be (see Figs. [3] 
and [TOD 

(ilnr 27rlnl0 o , . 

rwind~27r^^^ 1.45x10-3 s. (18) 

This timescale will correspond to the steep rising of energy of B^ around t = 0.1 s in 
Figs. [12] and [T3 After the steep growth, when the strength of Bfp becomes comparable 
to the poloidal comp onent, the growth ra te declines since B^^ grows by winding the 'weak' 



poloidal component (ITakiwaki et al. 20041 ). The final time of the simulation is determined 
when the Alfven speed reaches to the order of the speed of light. This is understood as 
follows: Alfven crossing time in the ^-direction at the innermost region becomes rA6/vA = 
10® X (7r/60)/c = 1.74 x 10~®s. Since the total timestep is several times of 10®, the final time 
is estimated to be several time, which is consistent with our results. The Alfven speed is 
estimated to be 

B B \ /10i2gcm-3V/' . ^ ^ 

VA = -y= ~ 2.82 X 10^ — — ^ cm s'^ (19) 
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which means that the final time is determined not by the time when the amphtude of 
reaches to lO^^G around the equatorial plane, but by the time when the amplitude of 
reaches to lO^^G at low density region, that is, around the polar region where the jet is 
launched. 

As for Br and Bq fields, from Fig. [T3l some instabilities seem to grow and saturate , 
which is similar to the behavior of MRI (IHawley fc Balbus 199ll : iBalbus fc Hawley 19981 ). 
At present, we consider that these instabilities are MRI modes with a wavelength of maximum 
growth mode unresolved. However, we can not conclude that these instabilities are really 
unresolved MRI modes. This is due to the reason as follows: The dispersion relation of the 
linear MRI modes is obtained analytically by assuming that the accretion disk is supported 
by rotation (that is, in Kepler motion). On the other hand, as shown in Figs. [3] and [101 the 
radial velocity is non-zero in the accretion disks in this study. Moreover, the radial flow speed 
is super-slow magnetosonic (the speed of slow magnetosonic wave is almost same with that 
of Alfven wave in this study: Fig. [TTl) . Since MRI is the instability of the slow- magnetosonic 
waves in a magnetized and differentially rotation plasma, the dispersion relation may be 
changed considerably for such a super-slow magnetosonic flow. However, there is no analytic 
solution for such a flow at present, so we use the dispersion relation of MRI for the discussion 
here. 



Ignoring entropy gradients, the condition for the ir istability of the slow-ni agnetosonic 
waves in a magnetized, differentially rotation plasma is (IBalbus &: Hawley 19911 ) 



dlnr 



+ (k-VA)' <0, 



(21) 



where k is the vector of the wave number. The wavelength of maximum growth of the linear 
instability is 



Ar 



27rt>A 



1.77 X 10^ 



lO^s 



B 



1013G 



lO^^g cm ^ 



1/2 



cm. 



(22) 



Since Aq is much smaller than the grid size of the inner most region (Ar = 3 x lO^cm), the 
linear MRI mode of maximum growth is not resolved in this study (note that the amplitudes 
of Br and Bq are of the order of lO^^G and much smaller than i?^; see Fig. fTT]) . However, MRI 
grows as long as Eq. fl2T]) holds. Since the value of the flrst term of Eq. fl2T]) at the inner most 
region is ~ —10^/ In 10 ~ —4.34 x 10'', it is conflrmed that inner most region is unstable for 
MRI mode with the wave lengt h longer than ~ lO^f-B/lO^^GldO^^g cm~^/p) Y^cm. The char- 
acteristic growing timescale is (IBalbus fc Hawley 19981 : lAkiyama et al. 20031 : iProga et al. 20031 ) 



tmri ~ 27r 



dhir 



-1/2 



6.74 X 10' 



(23) 
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which is seen in model 9 (Fig. [T3l) and model 8 (Fig. [T5l) . The saturation level of Br seems to 
be sli ghtly higher than that of Bg, which is similar to the results of the local simulations of 
MRI (jSano et al. 20041 ). Also, as shown in Fig. [121 the growth rate of B.^. and Bq components 
depends on the grid resolution, while the growth rate of B^ does not depend on the grid 
resolution. In fact, the growth rate becomes smaller for a coarse mesh case (Fig. [T6l(a)). 
This is similar to the picture that Br and Bg are amplified by MRI, while B^ is amplified 
by winding effect. The saturation level of Br and Bg seems not to be sensitive to the grid 
resolution. As stated above, in order to prove firmly that Br and Bg are amplified by MRI- 
like instability, dispersion relation of linear growing modes for super-slow magnetosonic flow 
has to be obtained analytically and the dispersion relation has to be reproduced by numerical 
simulations with finer grid resolution, which is outscope of this study. 

As shown in Fig. [TTl the estimated viscosity parameter becomes larger than 10"'^ at 
almost all region in r < 4 x 10^ cm. Since the angular velocity becomes larger than 10^ for 
r < 6 X 10^ cm (Fig. [TOl) . this viscosity becomes effective in a timescale of second, which 
is comparable to our simulations. Thus we think this viscosity drives high mass accretion 
rate in the models with magnetic fields. However, as stated in section 13.21 the angular 
momentum cannot be transfered to infinity along the radial direction. As shown in Fig. [TTl 
the flow is super- Alfvenic except for the inner most region. Thus angular momentum cannot 
be conveyed outward by Alfven wave in the radial direction. At present, we consider that 
the outflow (including the jet) in the polar direction in model 9 should bring the angular 
moment um from the inner most region (see Fig. [T] bottom). This picture is similar to 
CDAF (INarayan et al. 200ll ). We are planning to present further analysis of this feature in 
the forth-coming paper. 



Since the plasma beta can be of the order of unity at the inner most region (Fig. [141) 
and B^ dominates Br and Bg at the region (r < 7 — 8 x lO^cm), we can understand that 
the magnetic pressure from Br and Bg is much smaller than the radiation pressure (from 
photons, electrons, and positrons) and degenerate pressure. Thus, the jet cannot be launched 
by the effects of Br and Bg only. The winding-up effect is necessary to amplify i?^ field so 
that the magnetic pressure becomes comparable to the radiation and degenerate pressure. 
We saw that MRI-like instability seems to occur from the beginning of the simulations (see 
model 8 in Fig. [T^ . However, since the saturation level is not so high, the amplified energy 
of magnetic fields for Br and Bg cannot be seen well at first (see model 9 in Fig. [131 and model 
10 in Fig. [T5ll . This is because at the early phase the steady accretion disk is not formed 
and released gravitational energy is not so much (see Fig. [2]). In model 11 and model 12, 
the initial total energies of magnetic fields are so high that the effect of MRI-like instability 
cannot be seen (Fig. [151) . 
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As stated above, the axisymmetry of the system is assumed in this study So wind- 
ing up effect only amphfies m = mode of B^. Usually, it is pointed out that saturation 
level of the winding up effe cts becomes lower when three dimensional simulations are per- 
formed (IHawley et al. 19951 ). Thus 5^ may not be amplified as strong as lO^^G. Also, at the 
present study, the timescale for fields at the inner most region to be amplified to ~ lO^^G 
depends on the initial amplitude of the magnetic fields (see Figs. [13] and [15]) . However, if 
three dimensional simulations are performed, instabilities due to MRI(-like) modes do not 
decay due to Cowling's anti-dynamo theorem. Also, the MRI(-like) r nodes with wave vec - 



tors whose components are non-zero, which amplify Bfj, are included (iMasada et al. 2006[ ). 
If such MRI(-like) modes that amplify i?^ are included, the dependence of the dynamics 
on the initial amplitude of the magnetic fields may be diluted. We have to perform three 
dimensional calculations to see what happens in more realistic situations. 

As stated above, the timestep becomes so small when the Alfven crossing time becomes 
so small. Since this calculation is Newtonian, the speed of light is not included in the basic 
equations for macro physics (Eq. ([T])-Eq. (Q). In fact, we found that the Alfven speed 
becomes larger than the speed of light at some points by a factor of two or so at the final 
stage of simulations. Thus it will be one of the solutions to overcome this problem is to 
develop the special relativistic MHD code, in which the Alfven speed is, of course, solved to 
be smaller than the speed of light. 

In this study, we considered the ideal MHD without dissipation (it is, of course, numer- 
ical viscosity is inevitably included due to finite gridding effects). When resistive heating 
is efficient, considerable amount of energy in magnetic fields will be transfered to thermal 
energy by ohmic-like dissipation and reconnection, which will change the dynamics of collap- 
sars so much. The problem is, however, that the properties of resistivity of high density and 
high temperature matter with strong magnetic fields are highly uncertain (it is noted that 
artificial resistivity is included in Proga et al. (2003) in order to account of the dissipation 
in a controlled way instead of allowing numerical effects to dissipate magnetic fields in an 
uncontrolled manner (see also Stone and Pringle 2001)). 

Finally, we discuss the total explosion energy and bulk Lorentz factor of the jet. From 
Table [2] we can see that there seems to be a tendency that the mass of the jet becomes 
heavier when the initial amplitude of the magnetic fields are stronger (model 10, model 
11, and model 12), although this tendency is not monotonic (model 8 and model 9). This 
will reflect that the jet is launched earlier and the density is still higher for a case with 
stronger initial magnetic fields. There seems to be also tendency that the energy of the jet 
becomes also larger when the initial amplitude of the magnetic fields is set to be stronger, 
although this tendency is not also monotonic (model 8 and model 9). This tendency is not 
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so remarkable since the mass of the jet is heavier for a stronger magnetic field case and the 
mass should have some amount of kinetic, thermal, and magnetic field energies. As for the 
models 8 and 9, the information of initial condition might be lost considerably since it takes 
much time to launch jets. From these results, we can conclude that no GRB jet is realized 
even if strong magnetic fields is assumed. We consider that there will be a possibility that 
a GRB jet is realized if we can perform numerical simulations for much longer physical time 
(say, ~ 10 — 100s). In such long timescale simulations, the energy of fields should be 
much more than 10^° erg (see Figs. [T^ and [TH) . Also, the density in the jet may become lower 
along with time, because considerable mass will falls into the black hole along the polar axis. 
Thus the terminal bulk Lorentz factor may be enhanced at later phase. In order to achieve 
such a simulation, the special relativistic code will be helpful, as mentioned above. 



4.3. Prospect for Nucleosynthesis 

It is radioactive nuclei, ^^Ni and its daughter nuclei, ^^Co, that brighten the super- 
nova remnant and determine its bolometric luminosity. ^^Ni is considered to be synthesized 
through explosive nucleosynthesis because its half-life is very short (5.9 days). Thus it is 
natural to consider that explosive nucleosynthesis occurs in a hypernova that is accompanied 
by a GRB. However, it is not clearly known where the explosive nucleosynthesis occurs (e.g. 
Fryer et al. 2006a, Fryer et al. 2006b). 

Maeda et al. (2002) have done a numerical calculation of explosive nucleosynthesis 
launching a jet by depositing thermal and kinetic energy at the inner most region. They 
have shown that a mass of ^^Ni sufficient to explain the observation of hypernovae (~ O.5M0) 
can be synthesized around the jet region. In their calculation, all the explosion energy was 
deposited initially. Thus Nagataki et al. (2006) investigated the dependence of explosive 
nucleosynthesis on the energy deposition rate. They have shown that sufficient mass of 
^^Ni can be synthesized as long as all the explosion energy is deposited initially, while the 
synthesized mass of ^^Ni is insufficient if the explosion energy is deposited for 10 s (that is, 
the energy deposition rate is 10^^ erg s~^). This is because matter starts to move outward 
after the passage of the shock wave, and almost all of the matter moves away from the central 
engine before the injection of thermal energy (= 10^^ erg) is completed, so the amount of 
mass where complete burning to synthesize ^^Ni becomes little for such a long-duration 
explosion (see also Nagataki et al (2003)). 

On the other hand, it is pointed out the possibility that a substantial amount of ^^Ni is 
produced in the accreti on disk and a part of it is conyeyed outward by the viscosity-driven 



wind by some authors (IMacFadyen &: Woosley 19991 : iPruet et al. 20031 ). However, there is 
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much uncertainty how much ^^Ni is ejected from the accretion disk. This problem depends 
sensitively on the viscosity effects. Further investigation is required to estimate how much 
^^Ni is ejected. 

In this study, we have found that mass of ^^Ni in the accretion disk (see Fig. Hj) at the 
final stage of simulations is of the order of 10~^ Mq (M^°* in Table [3]). Thus if considerable 
fraction of the synthesized ^^Ni is ejected without falling into the black hole, there will 
be a possibility to supply sufficient amount of ^^Ni required to explain the luminosity of a 
hypernova. However, in the present study, the ejected mass of ^^Ni was found to be only 
10~^^ — 10~^M(7) [M^^ in Table [3]). This is because the entropy per baryon in the jet is so high 
(Fig. [T7]) that the light elements such as n, p, and He dominate in the jet. Thus we could not 
show that sufficient mass of ^^Ni was ejected from the accretion disk in the present study. 
There will be two possibilities to extract enough amount of ^^Ni from the accretion disk. 
One is that ^^Ni might be extracted efficiently from the accretion disk at later phase. Since 
typical temperature of the accretion disk will be lower when the mass o f the central black hol e 
becomes larger, the entropy per baryon in the jet will be decreased (INagataki et al. 20021 ). 
This feature should suggest that ^^Ni dominates in the jet component at the late phase. The 
other one is that some kinds of viscosities might work to convey the matter in the accretion 
disk outward efficiently. In MacFadyen and Woosley (1999), they included a- viscosity and 
showed that considerable amount can be conveyed. Since a viscosity is not included in 
this study, such a feature was not seen. However, when three dimensional simulations with 
higher resolution than in this study are done, much more modes of magnetic fields should 
be resolved, and some of them might be responsible for viscosities and work like a-viscosity. 

We discuss the possibility to synthesize heavy elements in coUapsars. We have shown 
that neutron-rich matter with high entropy per baryon is ejected along the jet axis (Figs. [17] 
and [TS]). This is because electron capture dominates positron capture in the accretion 
dis k (Fig. El). Thus there is a possib i lity that r-pro cess nucleosynthesis occurs in the 



Naeataki 2o"ool: iNagataki 200ll : [Nagataki fc Kohri 200ll : lwanaio et al. 2002 



Fujimoto et al. 20061 ). Moreover, we found that mass element with 



iet (INagataki et al. 1997 : 
Suzuki fc Nagataki 20oi 
high Ye (~ 0.522) appears around the polar axis near the black hole. This is because z/g cap- 
ture dominates Ue capture at the region. This is because flux of z/g from the accretion disk 
(note that z/g comes from electron capture) is sufficiently large to enhance Ye at the region. 
In su ch a high Yp region, there will be a possibihty that r/p— process nucleosynthesis oc- 
curs ( IWanajo 20061 ). We are planning to perform such a numerical simulations in the very 
near future. Also, there is a possibility that much neutrons are ejected from the jet since 
the entropy per baryon amounts to the order of 10^. If so, there may be a possibility that 
signals from the neutron decays may be observe d as the delayed bump of the light curve 
of the afterglow (IKulukarni 20051 ) or gamma-rays (IRazzaque fc Meszaros 2006al ). Also, GeV 
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emission may be observed by proton-neutron inelastic scattering (IMeszaros &: Rees 2000 



Rossi et al. 20061 : iRazzaque &: Meszaros 2006bl ) Finally, note that a careful treatment of 



neutrino transfer is important for the r- and r/p-process nucleosynthesis. As stated in sec- 
tion l4.1l the optically thin limit may break down at the high density region. When the matter 
becomes to be opaque to neutrinos, the neutrinos becomes to be trapped and degenerate. In 
the high density limit, the chemical equihbrium is achieved as /ie+/ip = fJ^n+fJ^ue where /le, fip, 
fj.n, and fii,^ are chemical potentials of electrons, protons, neutrons, and electron- type neu- 
trinos. When the chemical potential of electron-type neutrinos is not negligible, the electro n 
capture do not proceed further and electron fraction does not decrease so much (ISato 19751 ). 
which should be crucial to the r- and r/p-process nucleosynthesis. 



4.4. General Relativistic Effects 

In this section, we discuss general relativistic effects that we are planning to include in 
our magnetohydrodynamic simulation code. 

First, we still believe that effects of energy deposition due to weak interactions (espe- 
cially, neutrino anti-neutrino pair annihilation) can be a key process as the central engine 
of GRBs. In fact, the temperature of the accretion disk becomes higher especially for a 
Kerr black hole since the gravitational potential becomes deeper and much gravitational 
energy is released at the inner most re gion and the radius of the inner most stable orbi t 



becomes smaller for a Kerr black hole (jPopham et al. 19991 : iMacFadyen fc Woosley 19991 ). 



This effect will enhance the luminosity of the neutrinos from the accretion disk (Popham et 



al. 1999), and the neutrino pair an nihilation (lAsano fc Fukuyama 2000l : iMiller et al. 2003 



Kneller et al. 20061 : iGu et al. 20061 ). In the vicinity of the black hole, most of neutrinos and 
anti-neutrinos come shadows due to the bending effects of the neutrino geodesies. Since 
the shape and position of the black hole shadow depend on the physical parameters of the 
black hole (Bardeen 1973; Takahashi 2004, 2005), the effective area emitting the most of the 
neutrino flux just outside of the black hole shadows is also determined by the black hole 
parameters. These effects are included only when geodesies and disk structures are precisely 
calculated. 

Since according to our simulations the plasma beta (= pgas-|-radiation/pmag) in the 
polar region is lower than the plasma beta in other regions and the minimum values of the 
plasma beta is 0.193. This means that the polar region is magnetically dominated. In such 
regions where the magnetic energy is dominated such as a force-free field, the extraction of the 
rotational energy of the black hole is expected as Poynting flux (Blandford &; Znajek 1977) 
and the negative energy (Blandford & Payne 1982; Takahashi et al. 1990). Although these 
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effects have been c onfirmed numerically by many authors 



Komissarov 20051 : iHawley &: Krolic 20061 : iMcKinney 20061 ). there are many unsolved issues 



Koide 20031 : iMizuno et al. 2004bl : 



on the central engine of GRBs. Especially, the neutrino radiation from the negative energy 
fluid in the realistic general relativistic accretion flow is one of the unsolved issue. The 
past studies roughly estimate the extraction energy due to the Blandford-Znajek process to 
become non-neghgible amount to the explosion of the long GRSs (e.g. Lee et al. 2000a, 
2000b; Di Matteo et al. 2002). 



4.5. Variety of Initial Conditions 



The properties of a progenitor of GRBs are still unknown. There are much uncertainties 
on mass loss rate, final mass, angular momentum distribution, amplitude and configuration 
of magnetic fields prior to collapse. These properties will depend not only on initial mass and 
metall icity of progenitors but also on the presence (or absence) and properties o f companion 
stars jMaeder fc Mevnet 200l1 : IWooslev fc Heger 20061 : iMacFadven et al. 2005h . Thus it is 
necessary to perform numerical simulations with many models of progenitors to ensure the 
validity of the mechanism to launch a GRB jet from progenitors. In particular, we consider 
that the initial angular momentum distribution should be important since the mass accretion 
rate depends on it very sensitively. 



4.6. Prospects for Improvements of Numerical Scheme 

As shown in Eq. fl22|) . the critical wavelength of maximum growth of the linear MRI 
instability is quite short, so adaptive mesh refinement method (e.g. Norman 2005; Zhang 
and MacFadyen 2006; Morsony et al. 2006) is inevitably required to resolve the critical 
wavelength in a collapsar. Of course, much CPU time is required to perform three dimen- 
sional calculation with adaptive mesh refinement method. Moreover, when we try to include 
neutrino transfer code, three dimensions are required additionally for the momentum space. 
However, people will be able to realize such an expensive simulation in the near future with 
a help of rapidly growing power of super-computers. Of course, we are planning to devote 
ourselves to developing such a numerical code. 
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5. SUMMARY AND CONCLUSION 

We have performed two-dimensional magnetohydrodynamic simulations by the ZEUS- 
2D code to investigate the dynamics of a collapsar that generates a GRB jet, taking account 
of realistic equation of state (contribution of electrons, positrons, radiation, and ideal gas of 
nuclei), neutrino cooling and heating processes, magnetic fields, and gravitational force from 
the central black hole and self gravity. 

We have found that neutrino heating processes (neutrino and anti-neutrino pair annihi- 
lation, and fe and Pg captures on free nucleons) are not so efficient to launch a jet. We have 
found that a jet is launched by magnetic fields (in particular, fields that are amplified by 
the winding-up effect). However, the ratio of total energy relative to the rest mass energy in 
the jet at the final stage of simulations suggests that the bulk Lorentz factor of the jet will 
not reach to as high as several hundred, so we conclude the jet seen in this study will not be 
a GRB jet. We also found that the mass accretion rate seems to be enhanced in the models 
with magnetic fields. This might be because angular momentum is efficiently transfered by 
the viscosity due to the magnetic fields and multi-dimensional fiow. 

Since GRB jets are not obtained in this study, we consider that general relativistic ef- 
fects, by which the efficiency of energy deposition through weak interactions will be enhanced 
and rotation energy of the black hole will be transfered to matter through the magnetic fields, 
will be important to generate a GRB jet. Thus we are planning to develop a general rela- 
tivistic magnetohydrodynamics code in the very near future. Also, the accretion disk with 
magnetic fields may still play an important role to launch a GRB jet, and it may be seen 
if we can perform numerical simulations for much longer physical time (say, ~ 10 — 100 s). 
To reahze such a simulation, the special relativistic code will be helpful because the Alfven 
velocity is hmited to the speed of fight. 

We have shown that considerable amount of ^^Ni is synthesized in the accretion disk. 
Thus if some fraction of the synthesized ^^Ni is ejected without faUing into the black hole, 
there will be a possibihty for the accretion disk to supply sufficient amount of ^^Ni required 
to explain the luminosity of a hypernova. Also, we have shown that neutron-rich matter with 
high entropy per baryon is ejected along the rotation axis. This is because electron capture 
dominates positron capture. Moreover, we found that the electron fraction becomes larger 
than 0.5 around the polar axis near the black hole. This is because capture dominates 
Pg capture at the region. Thus there will be a possibility that r-process and r/p— process 
nucleosynthesis occur at these regions. Finally, much neutrons will be ejected from the jet, 
which suggests that signals from the neutron decays may be observed as the delayed bump 
of the light curve of the afterglow or gamma-rays. 
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Fig. 1. — Top panels: density contour with velocity fields for model at t = 2.1 s (left 
panel) and t = 2.2 s (right panel). Bottom panels: same with upper panels, but for model 
9. The color represents the density (g cm~^) in logarithmic scale (10^ — 10^^). In this figure, 
the central region of the progenitor (r < lO^cm) is shown. Vertical axis and horizontal axis 
represent rotation axis and equatorial plane, respectively. 
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Fig. 2. — Solid, short-dashed, dotted, dot-dashed hnes represent accreted total energy, kinetic 
energy, thermal energy, and electro-magnetic energy (xlO^° erg) into the black hole as a 
function of time, respectively. Long-dashed lines represent the mass accretion rate {Mq s~^) 
as a function of time. Left panel shows the case for model 0, while right panel represents 
the case for model 9. Note that in the case of model 0, the electro-magnetic energy is set to 
be zero. 
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Fig. 3. — Profiles of density, absolute value of radial velocity, angular frequency, temperature, 
density scale height, and specific angular momentum at the equatorial plane for model at 
t = 2.2 s. 
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Fig. 4. — Profiles of mass fraction for nuclear elements at the equatorial plane for model 
at t = 2.2 s. Dot-dashed, dotted, short-dashed, long-dashed, and solid lines represent mass 
fraction of n, p, ^He, ^^O, and ^^Ni, respectively. 
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Fig. 5. — Contour of neutrino cooling rate with velocity fields for model at t = 2.2 s. The 
color represents the emissivity of neutrinos (erg cm~^ s~^) in logarithmic scale (10^° — 10^^). 
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Fig. 6. — Cumulative energy (xlO^° erg) of emitted neutrinos for each process as a function 
of time. Dot-dashed, long-dashed, sohd, short-dashed, and dotted hues represent plasmon 
decay, electron-positron pair annihilation, positron capture, electron capture, and summation 
of all processes. Left panel shows the case for model 0, while right panel shows model 9. 
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Fig. 7. — Contour of neutrino heating rate with velocity fields for model at t = 2.2 s. The 
color represents the energy deposition rate (erg cm~^ s~^) in logarithmic scale (10^° — 10'^^). 
Left panel shows the energy deposition rate due to Ve and Ve captures on free nucleons, while 
right panel shows the energy deposition rate due to v and v pair annihilations (three fiavors 
are taken into account). 
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Fig. 8. — Neutrino luminosity (solid lines), energy deposition rate due to u and u pair 
annihilations (dotted lines), and energy deposition rate due to Uf, and captures on free 
nucleons (dashed hne) as a function of time. The unit is 10^° erg s~^. Left panel shows the 
case for model 0, while right panel shows the case for model 9. 
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Fig. 9. — Integrated deposited energy (xlO'''*' erg ) due to u and u pair annihilations (solid 
lines) and Ue and Ue captures on free nucleons (dashed lines) as a function of time. Left panel 
shows the case for model 0, while right panel shows the case for model 9. 
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Fig. 10. — Same with Fig.3, but for model 9. Note that the radial velocity profile is not 
similar to model at small radius. 
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Fig. 11. — Profile of amplitude of magnetic fields and estimated viscosity parameter a 
(xlO~^°) on the equatorial plane for model 9 at i = 2.2 s. Dotted, solid, and short-dashed 

lines represent the amplitude of B^, Bg, and B^, while long-dashed line represent the esti- 
mated viscosity parameter. Dot-dashed line represents the Alfven Mach number (xlO~^°) 
in the radial direction. 
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Fig. 12. — Profiles of angular momentum density (g cm~^ s~^) for model (long-dashed 
line) and model 9 (solid line) on the equatorial plane ai t — 2.2 s. Profiles of specific angular 
momentum (cm^ s~^) are also shown for model (dotted line) and model 9 (short-dashed 
hne) . 
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Fig. 13. — Evolution of total energy of magnetic fields (xlO^° erg ) for the case of model 9. 
Dotted, solid, and dashed lines represent energy in the form of Br, Bg, and B^, respectively. 
Note that 1-D simulation of the spherical collapse of the progenitor is done until t —0.1 s 
(see Section 2.2), so B^ is set to be until t — 0.1 s. 
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Fig. 14. — Contour of plasma beta (=Pgas+radiation/Pmag) with magnetic fields {Br and Be) for 
model 9 at t = 2.2 s. The color represents the plasma beta in logarithmic scale (10^^ — 10^). 
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Fig. 15. — Same as Fig. 13, but for model 8 (top- left panel), model 10 (top- right panel), 
model 11 (bottom-left panel), and model 12 (bottom right panel). 
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Fig. 16. — Left panel: Same as Fig. 13, but for 150(r) x5(6') grid points (a) and 150(r) x20(^) 
grid points (b). Right panel: Same as Fig. 13, but for 300(r)x60(^) grid points. In this 
simulation, neutrino anti-neutrino pair annihilation effect is not included to save CPU time. 
Simulation region is (lO^cm < r < 10^ cm; < ^ < 90). 
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Fig. 17. — Contour of entropy per baryon in units of kh a.t t = 2.2s for model 9. The range 
of the contour is from 1 to 10^. 



- 52 - 




Fig. 18. — Contour of electron fraction with velocity fields for for model at t = 2.2 s (left 
panel) and model 9 at t = 2.2 s (right panel). The color represents the electron fraction on 
a linear scale (0.1-0.540). 
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Model 


\\T/W\\ {%) 


Em/W\\ {%) 


Bo (G) 


Model 


1.3 








Model 8 


1.3 


1.1 X 10-s 


10« 


Model 9 


1.3 


1.1 X 10-6 


109 


Model 10 


1.3 


1.1 X 10"^ 


10^° 


Model 11 


1.3 


1.1 X 10-2 


10^1 


Model 12 


1.3 


1.1 


10^2 



Table 1: ||T/W^||: Initial ratio of the rotational energy to the gravitational energy. HE'm/W^lh 
Initial ratio of the magnetic energy. Bq is the strength of the magnetic field in the sphere 
(r < 3.6 X lO^cm). The digit in the name of each model represents the power index of Bq. 



Model 






Ei (erg) 


Mj (Mq) 


r/-i 




Model 8 {t 




= 5°) 


6.97E+47 


4.88E-6 


7.98E-2 


4.75E-2 


Model 8 {t 


h-- 


= 10°) 


7.40E+47 


5.39E-6 


7.67E-2 


5.36E-2 


Model 8 {i 


h-- 


= 15°) 


7.66E+47 


5.56E-6 


7.70E-2 


6.12E-2 


Model 9 {i 


h-- 


= 5°) 


2.96E+47 


2.42E-7 


6.81E-1 


4.71E-2 


Model 9 {t 


h-- 


= 10°) 


4.27E+47 


3.13E-7 


7.60E-1 


3.60E-1 


Model 9 (^j = 


= 15°) 


4.87E+47 


4.38E-7 


6.21E-1 


3.45E-2 


Model 10 


% 


= 5°) 


8.17E+45 


1.97E-8 


2.32E-1 


2.83E-1 


Model 10 




= 10°) 


9.08E+45 


2.12E-8 


2.39E-1 


2.72E-1 


Model 10 


[^J 


= 15°) 


1.32E+46 


3.65E-8 


2.03E-1 


2.13E-1 


Model 11 


[^J 


= 5°) 


1.05E+48 


1.16E-5 


5.05E-2 


2.74E-1 


Model 11 




= 10°) 


1.16E+48 


1.24E-5 


5.22E-2 


2.56E-1 


Model 11 


[^J 


= 15°) 


1.32E+48 


1.36E-5 


5.44E-2 


2.33E-1 


Model 12 


(^J 


= 5°) 


1.45E+49 


1.16E-4 


6.99E-2 


1.69E-1 


Model 12 


(^J 


= 10°) 


1.77E+49 


1.52E-4 


6.50E-2 


1.57E-1 


Mo (1(4 12 


^^.1 


= 15°) 


2.39E+19 


2.17E-1 


5.12E-2 


1.29E-1 



Table 2: (erg), Mj (Mq), F'^, and E^'^^/Ej are mass, total energy, terminal bulk Lorcntz 
factor, and ratio of the magnetic energy relative to total energy of the jet at the final stage 
of the simulations, respectively. 9j is the assumed opening angle of the jet. See section 3.3.1 
for details. 
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Model 




M^f (Me) 


Model 8 


6.04E-10 


3.79E-3 


Model 9 


7.52E-11 


5.92E-3 


Model 10 


1.09E-7 


3.95E-3 


Model 11 


1.40E-6 


1.61E-3 


Model 12 


5.62E-7 


1.12E-3 



Table 3: (Mq) and M^°* (Mq) represent the mass of ^^Ni in the regions where total 

energy (i.e. summation of kinetic energy, thermal energy, and gravitational energy) is positive 
and the total mass of ^^Ni in the whole simulated region at the final stage of the simulations, 
respectively. 



